#!/usr/bin/perl -w
##This script generated random bed/gtf file for an input bed/gtf file and from a selected region. 
##e.g. give an input of exon bed file, and the region of intron file, it will generate a random intron region for this exon bed file.
#cut -f1 knownGene-tss-UCSC062409-unique.1stExon.transcriptId.hg19.gtf | uniq -c | perl -ne 'chomp;if(/(chr\S+)/){$chr=$1;}print "$chr\n";' | less
## author: Yaping Liu  lyping1986@gmail.com 
## time: 2013-7-15

#Usage: perl randomIntervalFromSelectRegion.pl [option] input.bed/gtf/gff

use strict;
use Getopt::Long;
use File::Basename;

sub usage {
	print STDERR "\nUsage:\n";
    print STDERR "perl randomIntervalFromSelectRegion.pl [option] input.bed/gtf/gff output.bed/gtf/gff\n\n";
    print STDERR "This script generated random bed/gtf file for an input bed/gtf file and from a selected region.\n\n";
    print STDERR "  [Options]:\n\n";
	print STDERR "  --g FILE : A selected region file or whole genome interval.Should be sorted by coordinates and bed format. Required!\n\n";
	
	exit(1);
}
my $interval="";
GetOptions( "g=s" => \$interval,
			);
			
my $input=$ARGV[0];
my $output=$ARGV[1];
if(scalar(@ARGV)<2 || $interval eq ""){
	&usage();
}

open(FH,"<$input") or die "can not open file:$!";
#my $random_file = $input;
#$random_file =~ s/\.(\w+)$/.random.$1/;
open(OUT,">$output") or die "can not open file:$!";
my $pre_chr="";
my $tmp_interval="";
my $contig_range= "";
while(<FH>){
	chomp;
	my @splitin = split "\t";	
	
	if($pre_chr eq "" || $splitin[0] ne $pre_chr){
		$pre_chr=$splitin[0];
		if($tmp_interval ne ""){
			unlink($tmp_interval);
		}
		($contig_range,$tmp_interval) = &create_tmp_interval($pre_chr,$interval);
	}
	#print "$pre_chr\t$contig_range\t$tmp_interval\n";
	if($input=~/\.bed$/){
		my ($random_start,$random_end)=&creat_random_location($contig_range,$tmp_interval,$splitin[2]-$splitin[1]);
		$splitin[1]=$random_start;
		$splitin[2]=$random_end;
		if($splitin[2]-$splitin[1] > 0){ 
			my $line=join "\t",@splitin;
			print OUT "$line\n";
		}
	}
	elsif($input=~/\.gtf$/ or $input=~/\.gff$/){
		my ($random_start,$random_end)=&creat_random_location($contig_range,$tmp_interval,$splitin[4]-$splitin[3]+1);
		$splitin[3]=$random_start+1;
		$splitin[4]=$random_end;
		if($splitin[4]-$splitin[3] > 0){ 
			my $line=join "\t",@splitin;
			print OUT "$line\n";
		}
	}
	else{
		die "unrecoganized input file input format! Only allow bed/gtf/gff file format!\n\n";
	}
	
	
	
}
if($tmp_interval ne ""){
	unlink($tmp_interval);
}
close(FH);
close(OUT);

sub creat_random_location{
	my $contig_range = shift(@_);
	my $tmp_interval = shift(@_);
	my $block = shift(@_);
	

	my $random_number_start = int(rand($contig_range));
	my $start=$random_number_start;
	my $end=$start;

	open(GE,"<$tmp_interval") or die "cant open file $tmp_interval:$!\n";
	while(<GE>){
		chomp;
		my @splitin = split "\t";
		$start-=($splitin[2]-$splitin[1]);
		if($start<=0){
			$start+=$splitin[2]-$splitin[1];
			$start+=$splitin[1];
			$end=$start+$block;
			$end=$end<$splitin[2] ? $end: $splitin[2];
			last;
		}
	}
	close(GE);
	return ($start, $end);
}

sub create_tmp_interval{
	my $chr=shift @_;
	my $file= shift @_;
	my $tmp_interval=$file;
	$tmp_interval=~s/\.(\w+)$/.tmp_$chr.$1/;
	`grep "\^$chr\\b" $file > $tmp_interval`;
	my $contig_range= `awk '{SUM=SUM+(\$3-\$2)}END{print SUM}' $tmp_interval`;
	chomp($contig_range);
	return ($contig_range,$tmp_interval);
	
	
}